function  profitModel(m1,m2,sd1,sd2, r12, w1, w2, pw1, pw2,  N);

% pw1 = 10;    %price cap pay for ILEC for data
% pw2 = 10;    %price cap pay for ILEC for video

% Fixed Cost for ILEC
[CR1, CR2, CW1, CW2] = cost(N, N);

for Pi1=1:w1
    for Pi2=1:w2
        for Pc1=1:w1
            for Pc2=1:w2
                
                [p1, index1] = min(Pi1, Pc1);
                [p2, index2] = min(Pi2, Pc2);
                ProfitI(p1, p2) = 0;
                ProfitC(p1, p2) = 0;
                
                [surplus,n1, n2] = demandModel3D(m1,m2,sd1,sd2, r12, w1, w2, p1, p2,  N);
                CSurplus(p1, p2) = surplus;
                [cr1, cr2, cw1, cw2] = cost(n1, n2);
                
                %calculate the profit for product 1
                ProfitI\(p1, p2) = ProfitI(p1,p2) - N*(CR1 + CW1);
                
                if( index1 == 1 )
                    ProfitI(p1,p2) = ProfitI(p1,p2) + p1*n1;
                else           
                    if( pw1 < cw1 )
                        ProfitI(p1,p2) = ProfitI(p1,p2) + n1*pw1;
                        ProfitC(p1,p2) = ProfitC(p1,p2) - n1*pw1 - n1*cr1;
                    else
                        ProfitC(p1,p2) = ProfitC(p1,p2) - n1*(cw1+cr1);                
                    end
                        
                    ProfitC(p1,p2) = ProfitC(p1,p2) + p1*n1;
                end
                
                %calculate the profit for product 2
                ProfitI(p1, p2) = ProfitI(p1,p2) - N*(CR2 + CW2);
                
                if( index2 == 1 )
                    ProfitI(p1,p2) = ProfitI(p1,p2) + p2*n2;
                else                    
                    if( pw2 < cw2 )
                        ProfitI(p1,p2) = ProfitI(p1,p2) + n2*pw2;
                        ProfitC(p1,p2) = ProfitC(p1,p2) - n2*pw2 - n2*cr2;
                    else
                        ProfitC(p1,p2) = ProfitC(p1,p2) - n2*(cw2 + cr2);
                    end

                    ProfitC(p1,p2) = ProfitC(p1,p2) + p2*n2;
                end
                
            end
        end
    end
end

figure('NAME', 'CLEC Profit');
x=1:w1;
y=1:w2;
mesh(y,x,ProfitC)

figure('NAME', 'ILEC Profit');
x=1:w1;
y=1:w2;
mesh(y,x,ProfitI)


 figure('NAME', 'Consumer Surplus');
 mesh(y,x,CSurplus)
function [val, index] = min(p1, p2);
if(p1 < p2)
    val = p1;
    index = 1;
else
    val = p2;
    index = 2;
end